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Abstract. We present a general theoretical formulation to describe the interlayer interaction 
in incommensurate bilayer systems with arbitrary crystal structures. By starting from the 
tight-binding model with the distance-dependent transfer integral, we show that the interlayer 
coupling, which is highly complex in the real space, can be simply written in terms of generalized 
Umklapp process in the reciprocal space. The formulation is useful to describe the interaction in 
the two-dimensional interface of different materials with arbitrary lattice structures and relative 
orientations. We apply the method to the incommensurate bilayer graphene with a large rotation 
angle, which cannot be treated as a long-range moire superlattice, and obtain the quasi band 
structure and density of states within the first-order approximation. 
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1. INTRODUCTION 

Recently there has been extensive research efforts on atomically-thin nanomaterials, including 
graphene, hexagonal boron nitride(hBN), phosphorene and the metal transition dichalcogenides. 
The hybrid systems composed of different kind of atomic layers have also attracted much 
attention, and there the interaction at the interface between the different atomic layers often 
plays an important role in the physical property. In such a composite system, the lattice periods 
of the individual atomic layers are generally incommensurate due to the difference in the crystal 
structure and also due to misorientation between the adjacent layers. A well-known example of 
irregularly stacked multilayer system is the twisted bilayer graphene, in which two graphene 
layers are rotationally stacked at an arbitrary angle [Tj. When the rotation angle is small, 
in particular, the system exhibits a moire interference pattern of which period can be much 
greater than the atomic scale, and such a long-period modulation is known to strongly influence 
the low-energy electronic motion [HEllIlElEllDElia EH HD HS HS fl4] . Graphene-hBN 
composite system has also been intensively studied as another example of incommensurate 
multilayer system, where the two layers share the same hexagonal lattice structure but with 
slightly-different lattice constants, leading to the long-period modulation even at zero rotation 
angle HD Eg HD HU US ES ED- The electronic structure in graphene-hBN system was 

theoretically studied |22ll23ll21ll2ai2Sll2Zll2Sll22ll3iIll3DI32ll33ll33ll351, and the recent 

transport measurements revealed remarkable effects such as the formation of mini-Dirac bands 
and the fractal subband structure in magnetic fields [361 HD Ha [20]. 

The previous theoretical works mainly targeted the honeycomb lattice to model twisted 
bilayer graphene and graphene-hBN system, and also particularly focus on the long-period moire 
modulation which arises when the crystal structures of two layers are slightly different. Then 


we may ask how to treat general bilayer systems where the lattice vectors of the adjacent 
layers are not close to each other. In this paper, we develop a theoretical formulation to 
describe the interlayer interaction effect in general bilayer systems with arbitrary choice of 
crystal structures and relative orientations. By starting from the tight-binding model with 
distance-dependent transfer integral, we show that the interlayer coupling can be simply written 
in terms of a generalized Umklapp process in the reciprocal space. We then apply the method 
to the incommensurate bilayer graphene with a large rotation angle (6 = 20°) which cannot 
be treated as a long-range moire superlattice, and obtain the quasi band structure and density 
of states within the first-order approximation. Finally, we apply the formulation to the moire 
superlattice where the two lattice structures are close, and demonstrate that the long-range 
effective theory is naturally derived. 

2. Interlayer Hamiltonian for general incommensurate atomic layers 

We consider a bilayer system composed of a pair of two-dimensional atomic layers having 
generally different crystal structures as shown in Fig. [I] We write the primitive lattice vectors 
as ai and a 2 for layer 1 and ai and a 2 for layer 2, which are all along in-plane (x-y) direction. 
The reciprocal lattice vectors are defined by Gj and G* for layer 1 and 2, respectively, so as 
to satisfy a,; • G j = a.; • G j = 2n5ij. The area of the unit cell is given by S = |ai x a 2 | and 
S = |ai x § 2 1 for layer 1 and 2, respectively. Without specifying any details of the model, we 
can easily show that an electronic state with a Bloch wave vector k on layer 1 and one with k 
on layer 2 are coupled only when 

k + G = k + G, (1) 

where G = miGi + m 2 G 2 and G = rhiGi + ?7i2G2 are reciprocal lattice vectors of layer 1 and 2, 
respectively. This is regarded as a generalized Umklapp process between arbitrary misoriented 
crystals, and it can be easily understood by considering the wave decomposition in the plain 
wave basis as follows. A Bloch state on layer 1 (say is expressed as a summation of e*( k+G ) 
over the reciprocal vectors G, and one on layer 2 (</>? ) is expressed as a summation of e*( k+G ) 
over G. Also, Hamiltonian of the total system consists of Fourier components of G and G. As 
a result, the matrix element (4>^ \H\cj)^) can be non-zero only under the condition Eq. (fTj). 

In the following, we actually calculate the matrix elements for generalized Umklapp processes 
using the tight-binding model. We assume that a unit cell in each layer contains several atomic 

orbitals, which are specified by X = A, B, - • • for layer 1 and X = A, B, ■■ ■ for layer 2. The 

lattice positions are given by 

Rx = mai + n 2 a 2 + t x (layer 1), 

Ry = ni&i + n 2 £i 2 + (layer 2), (2) 

where n* and hi are integers, and tx and Ty are the sublattice position inside the unit cell, 
which can have in-plane and out-of-plane components. When the atomic layers are completely 
planar and stacked with interlayer spacing d, for instance, we have Tx ■ e z = 0 for layer 1 and 
Ty • e z = d for layer 2, where e z is the unit vector in z-direction. 

Let us define |Rx) = </>x(r — Rj) as the atomic state of the sublattice X localized at Rx- 
The atomic orbital 4>x may be different depending on X. We assume the transfer integral from 
the site Rx to R^ is written as —Ty X (Ry — Rx), i.e., depending on the relative position 
Rx — and also on the sort of atomic orbitals of X and X. The interlayer Hamiltonian to 
couple the layer 1 and 2 is then written as 

U = — ^ Ty_y(Ry ~ Rx) | Rjf ) (^ A' | T ^.C.. 


(3) 



Figure 1. Bilayer system composed of a pair of two-dimensional atomic layers with generally 
different crystal structures. Each unit cell can be generally composed of different sublattices 
and/or atomic orbitals. 


When the superlattice period is huge, constructing Hamiltonian in the real space bases becomes 
hard because it requires the relative inter-atom position for every single combination of the 
atomic sites of layer 1 and layer 2. On the other hand, the interlayer coupling is described in a 
simpler manner in the reciprocal space. We define the Bloch basis of each layer as 
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where k and k are the two-dimensional Bloch wave vectors parallel to the layer, N(N) is the 
number of unit cell of layer 1(2) in the total system area Stot- Although the layer 1 and layer 2 
are generally incommensurate, we assume the system has a large but finite area Stot = NS = NS 
to normalize the wave function. Then we can show that the matrix elements of U between Bloch 
bases can be written as, 

Ujt x { k,k) = (k,A|t/|k,A) 

= ~ ^xx(k + G) e * G TA+?6 Ta \_|_G,k+G’ ( 5 ) 
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which is non-zero only when the generalized Umklapp condition Eq. ([TJ is satisfied. Here t( q) 
is the in-plane Fourier transform of the transfer integral defined by 
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where z xx = (t x — tx ) • e z , and the integral in r is taken over the two-dimensional space of 5 tot . 
G and G run over all the reciprocal lattice vectors of layer 1 and 2, respectively. Since t xx (k) 
decays in large k, we only have a limited number of relevant components in the summation of 
Eq. ©. 

Eq. (J5]) is derived in a straightforward manner as follows. By inserting U in Eq. d3]) to the 
definition of U xx , we have 
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By applying the inverse Fourier transform of Eq. 0 , 
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to Eq. (ED, the second summation is transformed as 
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where we used in the last equation, 
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Using Eqs. ED and 0, we have 

Uxx( k. k) = - x E UxS + E 

G ri .Y 
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which is Eq. @. In the summation in R\- in the first line, we used the transformation 
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Figure 2. (a) Irregularly stacked bilayer graphene at rotation angle 8 = 20°. (b) Brillouin 

zones of the individual layers in the extended zone scheme. Symbols (filled circles, triangles and 
open circles) represent the positions of K + + G with several G’s. The different symbols indicate 
the different distances from the origin, (c) Corresponding positions of the symbols in (b) inside 
the first Brillouin zone of layer 2. 


3. Irregularly stacked honeycomb lattices 

We apply the general formulation obtained above to the irregularly stacked bilayer graphene 
with an arbitrary rotation angle. Here we consider a pair of hexagonal lattices as shown in Fig. 
[2|a), which are stacked with a relative rotation angle 8 and interlayer spacing d. We define the 
primitive lattice vectors of layer 1 as ai = o(l, 0) and a 2 = a(l/2, y/3/2) with the lattice constant 
a, and those of the layer 2 as a* = Ra^ where R is the rotation matrix by 8. Accordingly, the 
reciprocal lattice vectors of layer 1 (Gi, Go) and layer 2 are related by G, = RG*. The atomic 
positions are given by 

R.y = niai + n 2 a 2 + tx (layer 1), 

Ry = niai + n 2 a 2 + t^ (layer 2), (13) 

for X = A, B (layer 1) and X = A, B (layer 2), where 


ta = 0 , 

tb = -(ai + 2a 2 )/3, 
t a = de z + r 0 , 

+ r 0 - (^ + 2a 2 )/3. (14) 


Here we take the origin at an A site, and define To as the relative in-plane translation vector of 
the layer 2 to layer 1. 

To describe the electron’s motion, we adopt the single-orbital tight-binding model for p z 
atomic orbitals. Then Ty Y (R) does not depend on indexes X and X, and it is approximately 
written in terms of the Slater-Koster parametrization as, m 
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The typical parameters for graphene are a ~ 0.246nm, d « 0.335 nm, V pp7T « —2.7eV, 
V ppa « 0.48eV, and ro ~ 0.184a[ITj. Once the transfer integral T(R) between the atomic sites 
is given, one can calculate the in-plane Fourier transform f(q), and then obtain the interlayer 
Hamiltonian by Eq. ©• Since the transfer integral is isotropic along the in-plane direction, we 
can write t(q) = t(q ) with q = |q|. 

Fig. H a ) illustrates the lattice structure at rotation angle 6 = 20°. Fig. HL) shows the 
Brillouin zones of the two layers in the extended zone scheme, where the blue symbols (filled 
circles, triangles and open circles) represent the positions of k + G for some particular k (here 
chosen as the zone corner K + ) with several G’s. Figs. [21(c) plots the corresponding positions of 
these symbols inside the first Brillouin zone of layer 2. This indicates the wave numbers of layer 
2, k = k+G —G, which are coupled to k of layer 1 under the condition of Eq. ©• The amplitude 
of the coupling is given by t(k + G), and it solely depends on the distance to each symbol from 
the L-space origin in Fig. [2](b). In the present parameter choice, the amplitudes for filled circles, 
triangles and open circles are t{I\) ~ 110 meV, t(2I\) ~ 1.6 meV, and t(\/7K) ~ 0.062 meV 
respectively, where K = |K + | = 47r/(3a). The couplings for other fc-points are exponentially 
small and negligible. 

When the lattice vectors of the two layers are incommensurate (i.e., do not have a common 
multiple) as in this example, we cannot define the common Brillouin zone nor calculate the 
exact band structure, since the interlayer interaction connects infinite number of fc-points in the 
Brillouin zones of layer 1 and layer 2. But still, we can obtain an approximate band structure, 
considering only the first-order interlayer processes while neglecting multiple processes. Let 
us consider a particular wave vector k of layer 1, and take all k’s in layer 2 which are directly 
coupled to k as in Fig. He). By neglecting exponentially small matrix elements, we can construct 
a finite Hamiltonian matrix including only the bases of a single wave vector k of layer 1 and 
several k’s in layer 2. By diagonalizing the matrix, we obtain the energy eigenvalues £„k labeled 
by the index n. We then define the spectral function contributed from layer 1 as 

^i(k,e) = ~ £ ™k), (16) 

n 


where is the total wave amplitudes on layer 1 in the state |nk). The density of states 
contributed from layer 1 is expressed as 

Di(e) = ( 17 ) 

k 

where the summation is taken over the first Brillouin zone of layer 1. We perform the exactly 
same procedure for the layer 2 by considering a single k of layer 2 and all k’s in layer 1 coupled 
to k, and obtain the spectral function and the density of states for the layer 2. The total density 
of states of the system is given by D\ + D -2 . 

Fig.[3]plots the total density of state D(e) calculated for the incommensurate bilayer graphene 
with 6 = 20°. The red dashed curve is the density of states of decoupled bilayer graphene 
(i.e., twice of monolayer’s). We actually see additional peak structures due to the interlayer 
interaction, and these features are consistent with the exact tight-binding calculation for a 
commensurate bilayer graphene at a similar rotation angle m- Fig. [4] illustrates the Fermi 
surface reconstruction at several different energies. The right panel in each row shows the 
layer l’s spectral function Hi(k, e) in presence of the interlayer coupling. The left panel shows 
the Fermi surface in absence of the interlayer coupling, where the black curves represent the 
equienergy lines for the layer l’s energy dispersion £i(k), and the pink curves are for the layer 
2’s dispersion with &:-space shift, £ 2 (k + G). Thickness of the pink curves indicate the absolute 
value of the interlayer coupling t(k + G). 



Figure 3. Electronic density of states in the incommensurate bilayer graphene with 8 = 20°, 
calculated in the first-order approximation (see the text). The red dashed curve plots the density 
of states of decoupled bilayer graphene. Vertical arrows indicate the energies considered in Fig. 

m 


Fig. [H(a) shows a typical case (e = 2.5 eV) where we observe small band anticrossing at 
the intersection of the Fermi surfaces of layer 1 and layer 2. Fig. 0fb) (e = —2.85 eV) is 
for the energy at which the density of states exhibits a dip [Fig. [3] . There the interlayer 
coupling is relatively strong, and indeed we see that the original triangular Fermi pockets of the 
individual layers are strongly mixed and reconstructed into a large single Fermi surface. Note 
that the coupling strength depends not only on t(k + G), but also the relative phase factors 
between different sublattices. At even lower energy e = —5 eV [Fig. [2(b)] beyond the van-Hove 
singularity, the layer 1 and the layer 2 have almost identical Fermi surface surrounding T point 
in absence of the interlayer coupling, and they are hybridized into a pair of circles with different 
radii corresponding to the bonding and anti-bonding states. This feature is also observed in the 
density of states [Fig. [3] as a large split of the band bottom. Such a splitting does not occur 
in the positive energy region, because there A and B sublattices in the same layer have the 
opposite signs in the wave amplitude, so that the interlayer mixing vanishes due to the phase 
cancellation. 

In the above approximation, we neglect the second order processes such that k points of layer 
2 (linked from initial k of layer 1) are further coupled to other k 7 of layer 1. We can include such 
higher-order processes up to any desired order as follows. To include the second order process 
in calculating Ai(k, e), for example, we take all k of layer 2 which are coupled to k of layer 1 
and also take all k 7 of layer 1 which are coupled to k, and then construct a Hamiltonian matrix 
with a larger dimension. From the obtained eigenstates, we calculate the spectral function using 
Eq. (fTHl) . but then g ri 1 ^ should be the total wave amplitudes from layer 1 at the original k (0-th 
order), without including k 7 . 

4. Long-period moire superlattice 

When the primitive lattice vectors of layer 1 and those of layer 2 are slightly different, the 
interference of two lattice structures give rise to a long-period moire pattern, and then we can 
use the long-range effective theory to describe the interlayer interaction. in Gam u ei 0 in i 
uni on im iia nn Here we show that the long-range effective theory can be naturally derived 
from the present general formulation, just by assuming that the two lattice structures are close 
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Figure 4. Fermi surface reconstruction in the incommensurate bilayer graphene with 8 = 20°, 
at energies of (a) 2.5eV, (b) —2.85eV and (c) —5eV. The right panel in each row shows the layer 
l’s spectral function x4i(k, e) in presence of the interlayer coupling. The left panel shows the 
Fermi surface in absence of the interlayer coupling, where the black curve is for the layer l’s 
energy dispersion £i(k), and the pink curve is for the layer 2’s dispersion with fc-space shift, 
£o(k + G). Thickness of the pink curve indicates the absolute value of the interlayer hopping 
t(k + G). 
















Figure 5. Similar plots to Fig. [ 2 ] for 6 = 5°. 


to each other. 

We define a linear transformation with a matrix A that relates the primitive lattice vectors 
of layer 1 and 2 as 


ai = Aa { . (18) 

Correspondingly, the reciprocal lattice vectors become 

Gj = (A*)- 1 G i: (19) 

to satisfy a* • G j = aj • G j = 2%S l j. When the lattice structures of the two layers are similar, 
the matrix A is close to the unit matrix. The reciprocal lattice vectors of the moire superlattice 
is given by small difference between Gj and Gj as 

Gf = G i -G i = [l-(^ t )- 1 ]G i . (20) 

When the two layers are identical and rotationally stacked with a small angle, for example, the 
matrix A is given by a rotation R. Noting (Rt)- 1 = R , we have 

Gf = (l -R)Gi. (21) 

When layer 1 and layer 2 have different lattice constant as in the graphene-hBN bilayer, the 
matrix A is given by the combination of the isotropic expansion and the rotation as aR, where 
a is the lattice constant ratio. Since [(ai?)^ 1 = a~ 1 R, we have 

Gf = [1 - a~ 1 R]G i . (22) 

When the specific form of the matrix A is given, we can immediately derive the interlayer matrix 
elements for the long wavelength components using Eq. © as, 

Ux X ( k + miGi 1 + m 2 G^, k) = —t xx ( k + miGi + m. 2 G 2 ) 

x e -i(miGi+m 2 G2)-Tx+i(miGi+m2G2)-r i ^3) 


where mi and m 2 are integers. 



























































For example, let us derive the interlayer Hamiltonian of the irregularly stacked bilayer 
graphene with a small rotation angle. Since the low-energy spectrum of graphene is dominated 
by the electronic states around the Brillouin zone corners K and K' , we consider the matrix 
elements for initial and final £;-vectors near those points. The I\ and K' points are located at 
= — £(2Gi + G 2 )/2> for layer 1 and = — £(2Gi + G2)/3 for layer 2, where £ = ±1 are the 
valley indexes for K and K' , respectively. When we start from k = K + , for example, a typical 
scattering process is illustrated in Figs. 0(b) and (c) in a similar manner to Fig. 0 There the 
electron at K + in layer 1 is coupled to K + + m\G\ + m 2 G 2 in layer 2 with absolute amplitude 
t(K + +miGi+m 2 + G 2 ). As already argued in the previous section, the coupling amplitudes for 
filled circles, triangles and open circles in Figs. 0(b) and (c) are t(K) ~ 110 meV, t(2K) ~ 1.6 
meV, and t{\flK) ~ 0.062 meV respectively, and the couplings to other fc-points are negligibly 
small. The matrix element changes when the initial vector k is shifted from K+, but we neglect 
such a dependence assuming k is close to K + . As a result, we obtain the interlayer Hamiltonian 
of near from Eq. (1231) as 
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where r is the in-plane position, co = exp(27ri/3) and To is set to 0. 
The total Hamiltonian is written in the basis of {A, B, A, B} as 

U®-( H 1 

Aff - l U H 2 ’ 


(24) 


(25) 


where H\ and H 2 are the intralayer Hamiltonian of layer 1 and 2, respectively, defined by 


H 1 « -hv {k - K ? ) • (€a x ,a y ), 

H 2 « -hv[R ~\k - K^)] • (£<t x , a y ), (26) 

with Pauli matrices a x and a y , and the graphene’s band velocity v. If we only take the terms 
with t(K ) in the matrix U, the expression becomes consistent with the previous formulation 
for graphene-graphene bilayer. GQ El El DU Note that the Hamiltonian matrix depends on the 
actual choice of K points out of the equivalent Brillouin zone corners; + miGi + m 2 G 2 and 
Kjj + fhiGi + fh 2 G 2 . The different choice of ra; and rhi adds extra phase factors to the Bloch 
bases depending on sublattices, while the resulting Hamiltonian matrix are connected to the 
original by an unitary transformation. 

While we neglected the interlayer translation vector to in deriving Eq. (1241) . the matrix 
elements of U actually depend on t 0 according to Eq. (1231) . where the term with e*( miG i + m 2 G 2 )-r 
is accompanied by an additional phase factor g i (™i G i+ m 2 G 2 )-To i This extra term, however, can 
be incorporated into a shift of the space origin as 

6 i(miGf+m2G“)-r e «(miGi+m2G2)-ro _ gipraiG^+rr^G^Mr+ro) ^7) 

where ro = (A — 1)to- This reflects the fact that relative sliding between two layers leads to 
a shift of the moire interference pattern in the real space. The only exception is when the 






two layers shares the same lattice vectors Gj = Gj, where Gl' 1 vanishes so that To cannot 
be eliminated by shifting the origin. In this case, the energy band actually becomes different 
depending on the sliding vector To- For the graphene bilayer case, in particular, the expression 
of U with Gj = Gj becomes equivalent to that of regularly-stacked graphene bilayer with a 
interlayer sliding. [38] H9j [JOJ 

5. Conclusion 

We theoretically studied the interlayer interaction in general incommensurate bilayer systems 
with arbitrary crystal structures. Using the generic tight-binding formulation, we demonstrate 
that the interlayer coupling in the reciprocal space is simply expressed in terms of a generalized 
Umklapp process. We applied the formula to the incommensurate honeycomb lattice bilayer with 
a large rotation angle, which cannot be treated as a long-range moire superlattice, and actually 
obtain the quasi band structure and density of states within the first-order approximation. 
Finally, we apply the formulation to the moire superlattice where the two lattice structures are 
close, and derive the long-range effective theory with a straightforward calculation. 
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